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Abstract 

This paper presents a combined numerical-theoretical study of the macroscopic behavior and local field dis- 
tributions in a special class of two-dimensional periodic composites with viscoplastic phases. The emphasis 
is on strongly nonlinear materials containing pores or rigid inclusions. Full-field numerical simulations are 
carried out using a Fast-Fourier Transform algorithm [H. Moulinec, P. Suquet, C. R. Acad. Sci. Paris II 
318, 1417 (1994)], while the theoretical results are obtained by means of the 'second-order' nonlinear ho- 
mogenization method [P. Ponte Castaheda, J. Mech. Phys. Solids 50, 737 (2002)]. The effect of nonlinearity 
and inclusion concentration is investigated in the context of power-law (with strain-rate sensitivity m) be- 
havior for the matrix phase under in-plane shear loadings. Overall, the 'second-order' estimates are found 
to be in good agreement with the numerical simulations, with the best agreement for the rigidly reinforced 
materials. For the porous systems, as the nonlinearity increases (m decreases), the strain field is found 
to localize along shear bands passing through the voids (the strain fiuctuations becoming unbounded) and 
the effective stress exhibits a singular behavior in the dilute limit. More specifically, for small porosities 
and fixed nonlinearity m > 0, the effective stress decreases linearly with increasing porosity. However, for 
ideally plastic behavior (m = 0), the dependence on porosity becomes non-analytic. On the other hand, for 
rigidly-reinforced composites, the strain field adopts a tile pattern with bounded strain fiuctuations, and no 
singular behavior is observed (to leading order) in the dilute limit. 
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1. Introduction 

The overall (visco-) plastic fiow of composites exhibiting strong constitutive nonlinearities and large 
heterogeneity contrasts is particularly difficult to model due to the tendency of the strain (-rate) field to 
localize in shear bands running through the composite along 'minimal paths' (Drucker, 1966; Duxbury et al., 
2006). As full-field numerical simulations of such material systems demand significant computational power, 
especially when the microstructure is random, nonlinear effective-medium theories, also known as homoge- 
nization methods, are particularly attractive modeling tools. Most of the available theories rely on the use 
of a 'linear comparison medium' of one sort or another (see, for instance, reviews by Ponte Castaheda and 
Suquet (1998); Willis (2000)). Two approximations are typically involved in this class of theories. The first 
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one consists in linearizing the nonlinear behavior of the individual constituents following some appropriate 
linearization procedure, and the second one consists in estimating the overall response of the resulting linear 
comparison medium. For example, Hill (1965) proposed an incremental procedure making use of the tangent 
moduli of the phases evaluated at the average field in the phases and of the self-consistent linear homog- 
enization method to estimate the macroscopic response of elastoplastic poly crystals, whereas Hutchinson 
(1976) proposed the use of the secant modulus, also evaluated at the phase averages, together with the 
self- consistent method to estimate the response of viscoplastic polycrystals. The linearization schemes in 
both of these approaches were physically plausible, but ad hoc. Ponte Castaheda (1992) introduced new 
variational principles where the trial field was the moduli of a suitably chosen 'linear comparison composite.' 
By choosing constant moduli in the phases, this procedure leads to improved analytical estimates (Ponte 
Castaheda, 1991) for nonlinear composites, and as shown by Suquet (1995), the resulting linearization can be 
interpreted in terms of the secant moduli of the phases evaluated at the second moments of the fields in the 
phases. A more general procedure, making use of a 'linear thermoelastic composite' and having the property 
of delivering estimates that are exact to second order in the contrast, was developed by Ponte Castaheda 
(2002a). In addition, a related procedure, making use of a similar linearization, but exploiting instead the 
properties of an assumed Gaussian distributions for the fields, was advanced by Pellegrini (2001a). 

Recent comparisons with numerical simulations for two-phase materials (Idiart et al., 2006b; Rekik et 
al., 2007; Idiart and Ponte Castaheda, 2007b) and polycrystals (Lebensohn et al., 2004a, b, 2007) indicate 
that, within that class of 'linear comparison' theories, the so-called 'second-order' theory of Ponte Castaheda 
(2002a) delivers the most accurate predictions for the macroscopic response as well as for low-order statistics 
of the local fields. This theory has been found to capture some of the essential features associated with 
strain localization. For instance, it correctly predicts the anisotropic and unbounded character of the strain 
fluctuations in ideally plastic, particle-reinforced random composites (Idiart et al., 2006b; Idiart and Ponte 
Castaheda, 2007b), and the non-analytic dependence on porosity of the yield stress of random porous media 
with vanishingly small porosities (Ponte Castaheda, 2002a; Idiart and Ponte Castaheda, 2007b). 

The 'second-order' theory makes use of a linear comparison composite (LCC) with anisotropic phases, 
where the direction and strength of anisotropy is determined by the means and variances of the mechanical 
fields in each phase, and the constitutive nonlinearity. This work is aimed at understanding how the 'second- 
order' theory — and 'linear comparison' theories more generally — copes with strain localization in nonlinear 
composites by exploiting the constitutive anisotropy in the underlying LCC. To that end, this work focuses 
on a special class of tif o-dimensional periodic composites for which numerical simulations are simple enough 
so that a thorough and systematic comparison is feasible. The paper builds on the recent works of Willot 
et al. (2008a, b), which provided similar analyses and comparisons for strongly anisotropic linear elastic 
composites. As will be seen below, the anisotropy of the LCC in the context of the second-order estimates 
for the nonlinear composites in some sense plays the role of the nonlinearity, with strong nonlinearity leading 
to strong anisotropy in the LCC. The full-field numerical simulations are carried out using the Fast-Fourier 
Transform (FFT) algorithm of Moulinec and Suquet (1994) in the form recently proposed by Willot and 
Pellegrini (2008). We focus on porous materials and rigidly-reinforced composites, which are the cases of 
extreme heterogeneity contrast, and therefore good case studies. In addition, particular attention will be 
given to the dilute limits, which serve to illustrate the significant differences in behavior between these two 
extreme cases. 

2. Preliminaries on periodic composites 

2.1. Local and effective behavior 

The materials considered are made up of a continuous matrix phase (r = 1) containing cylindrical 
inclusions (r = 2) with circular cross section, aligned with the xa-axis and periodically distributed in the 
transverse plane as shown in Fig. 1. The inclusion phase will be taken as either vacuous or rigid. 

The constitutive behavior of the matrix phase is taken to be characterized by an isotropic^ incompressible 
strain potential w^^\ such that the stress and strain tensors are related by 

'^=^(^)' ^^^\^) = <t>{^e), (1) 
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Figure 1: Microstructure of the periodic material with reference axes and unit cell, together with a schematic representation 
of 'simple shear' (SS) and 'pure shear' (PS) macroscopic loadings. 



where the von Mises equivalent strain is defined in terms of the deviatoric part of the strain tensor by 
£e = {2/3)ed • and tr(6) = 0. This constitutive relation can be used within the context of the 
deformation theory of plasticity, where s and cr represent the infinitesimal strain and stress, respectively. 
Relation (1) applies equally well to viscoplastic materials, in which case e and a represent the Eulerian 
strain rate and Cauchy stress, respectively. A relation completely equivalent to (3) results from a dual 
formulation, which makes use of a stress potential u^^\ such that 

e=^{a), ti«{<r)=^(ae), (2) 

OCT 

where the von Mises equivalent stress is given in terms of the deviatoric stress tensor by cTg = \/ {3/2)crd - (Td- 
Let (•) denote the volume average over the composite. The effective behavior of the composite is defined 
as the relation between the average stress a = (cr) and the average strain s = (s) . When the ratio between 
the characteristic length scale of the unit cell to that of the specimen tends to zero, there exists (under 
certain technical assumptions) a homogenized behavior, which can be characterized by 

_ dw _ du 

OS da 

where w and u are the effective strain and stress potentials^ defined by 

w{^) = min_ (i(;(x, s)), 2(a) = min_ {u{^^cr))^ (4) 
6e/C(6) cres{cr) 

with lC{e) and 5(a) denoting the set of 'kinematically admissible' strain and 'statically admissible' fields 
(see Ponte Castafieda and Suquet (1998) for details). Thus, the problem of finding the effective behavior of 
the composite reduces to that of computing the effective potential w or u. 

For definiteness, it is assumed that the behavior of the matrix follows a power law of the form 

(\ l+m / \ 1+n 

'A , ^K) = f^f^) , (5) 

where ctq and Sq are the flow stress and the reference strain, respectively, n = 1/m denotes the nonlinearity, 
and m denotes the so-called strain-rate sensitivity, and is such that < m < 1. This model is particularly 
appropriate to study a wide range of material behaviors, such as the time-independent plastic deformation 
of metals as well as their time- dependent viscous deformation (i.e., high temperature creep). The limiting 
values m = 1 and m = correspond to linear and rigid-ideally plastic behaviors, respectively. 

We restrict the analysis to isochoric, plane-strain macroscopic deformations (or purely deviatoric stress). 
Because of the intrinsic anisotropy of its microstructure, the effective behavior of the composite material 
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will depend on the loading angles and ?/5 that are defined by the orientation of the principal axes of the 
deformation and loading, respectively, with the symmetry axes of the microstructure. We will focus on 
the two "extreme" cases of loading orientation: {i) 'pure shear' (PS) corresponding to the case where the 
principal axes of 6 and a coincide with the Cartesian axes introduced in Fig. 1, and (ii) 'simple shear' (SS) 
corresponding to the case where the principal axes of these tensors form 45° with respect to the Cartesian 
axes, as depicted in Fig. 1. It then follows from the homogeneity of the potentials (5) that, for these specific 
macroscopic loadings, the eff"ective potentials (3) can be written as 

~ /_ \ l+m ~ / — \ 1+n 

^ ^ l + m V^oy ^ ^ 1 + n Vcroy ^ ^ 



where Se = (2/v^)a/5i2 + (l/4)(6ii — 622)^ and = V^a/^i2 + (l/4)(o'ii —^22)^ are the macroscopic 
equivalent strain and stress, respectively, and 5o denotes the effective flow stress of the material, which 
depends on the specific macroscopic loading (SS or PS), the strain-rate sensitivity m and the inclusion 
volume fraction /. 

2. 2. Statistics of the local fields 

In addition to the effective behavior, homogenization methods can deliver useful information about the 
field distributions in the form of low-order statistics such as the first and second moments of the fields in 
each constituent phase (Idiart and Ponte Castaheda, 2007a). For the macroscopic loadings considered here, 
it is natural to express the deviatoric parts of the local fields as 

v^S 1 

= ^ (^ps eps + ^ss ess) , ^ ^ ^^^^ ~^ ^^^^ ' 

where the unit second-order tensors eps and ess are given in terms of the Cartesian axes of Fig. 1 by 

ess = ei (g) e2 + e2 (8) ei, eps = ei (g) ei - e2 (8) 62- (8) 

The prefactors in (7) are such that = + 6pg and similarly for the stress. We also introduce the 
fourth-order projection tensors (see Willot et al. (2008a)): 

= ^ess ^ ess, = ^eps eps- (9) 

Let (•)^^) denote the volume average over phase r. We are interested in the first moments (phase averages), 
denoted by 6^^^ = (s)^^^ and a^^^ = {ct)^^\ and the second moments (s s)*^^^ and {a (g) cr)^^^ of the field 
distributions. A measure of the field fluctuations within each phase is given by the phase covariance tensors 

= {e^ s)^^) - 6^^) (10) 

and similarly for the stress. For power-law composites like the ones considered here, it can be shown that the 
local strain and stress fields are homogeneous functions of degree 1 in Se and ^e, respectively. In addition, 
it follows from symmetry considerations that, for 'aligned' loadings in the sense described above, the phase 
averages of the local fields should be co- axial with the macroscopic averages. Then, the deviatoric part of 
the phase averages can be written as 



_(r) _(r) 
_(r) — — (r) — /-i -i \ 

^d-^^d, cry = —crd, (11) 



where the ratios e^e^ jse and cr^e^ /cje depend on the strain-rate sensitivity, the porosity, and the loading 
direction (PS or SS). In turn, it also follows from symmetry considerations that the corresponding second 
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moments should be 'aligned' with the macroscopic averages, in the sense that one of their eigentensors is co- 
axial with e and a. Then, the standard deviations {SD^'^\-) = \/ ((O^)*^^^ — ((')*^^^)^) ^^^^ components 
defined in (7) are obtained by projecting the covariance tensors (10) onto the tensors (9): 

5DW(ess) = y^ESS::4'-\ 5D«(ass) = ^^ESS ::C«, (12) 

and similarly for the PS components. For convenience, when discussing the results, we will refer to the 
SS and PS components of the strain/stress tensors as 'parallel' and 'perpendicular' to the macroscopic 
strain/stress tensors in the case of SS loadings, and conversely in the case of PS loadings. 

2.3. Duality between porous materials and rigidly-reinforced composites 

Table 1: Duality relations between power-law, two-dimensional materials with pores and rigid inclusions. The tensor Q 
represents a 7r/4 counterclockwise rotation in the plane. 
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The methods to be described in Sections 2.4 and 2.5 are general and can be directly applied to porous 
materials and rigidly-reinforced composites. In this work, however, we exploit the duality relations of 
Francfort and Suquet (2001) for incompressible, nonlinear, two-dimensional elasticity (see Proposition 5.3 
in that Ref.), which allow us to obtain results for rigidly-reinforced materials directly from corresponding 
results for porous materials. In the case of a power-law matrix, results for a rigidly-reinforced composite with 
a strain-rate sensitivity m can be obtained from corresponding results for a porous material with strain-rate 
sensitivity l/m using the relations given in Table 1. Note that we are therefore required to solve for porous 
materials with strain-rate sensitivities larger than 1. Thus, the methods need to be implemented for porous 
materials only. 

2.4' Second- order variational method 

The so-called 'second-order' method is a fairly general method for estimating the effective potentials 
of nonlinear composites, introduced by Ponte Castaheda (2002a), which delivers estimates that are exact 
to second order in the heterogeneity contrast. The central idea of this method is to introduce a linear 
comparison composite (LCC), with the same microstructure as the nonlinear composite, whose constituent 
phases are identified with appropriate linearizations of the given nonlinear phases, determined through a 
variational procedure. This allows the use of the many different methods already available to estimate the 
effective potentials of linear composites to generate corresponding estimates for the effective potentials of 
nonlinear composites. Estimates for the effective behavior then follow from relations (3), and estimates for 
the field statistics follow from similar identities (Pellegrini, 2001b; Idiart and Ponte Castaheda, 2007a). 

The 'second-order' method makes use of a LCC with local elasticity tensors given by generalized secant 
moduli that are intermediate between the standard secant and tangent moduli of the nonlinear phases, 
employed in earlier methods. In general, these elasticity tensors are anisotropic, even if the nonlinear phases 
are isotropic. The direction and strength of anisotropy are determined by the means and variances of 
the mechanical fields in the LCC, and the constitutive nonlinearity. For the particular class of material 
systems and loading conditions considered in this work, the matrix phase in the LCC is incompressible and 
characterized by two different shear moduli A and /i, one for each of the shear deformation modes defined 
in (8). The degree of anisotropy is thus given by the ratio k = \/ (i between the two moduli, with the 
extreme values k = 1 and A: = corresponding to isotropic and strongly anisotropic phases, respectively. 
Linear materials of this type have been studied in detail in Willot et al. (2008a). The linear estimates of the 
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Hashin-Shtrikman type derived in that paper were found to be very accurate, relative to FFT simulations, 
for all values of k and inclusion concentrations in the range < / < 0.3. Those estimates are used here to 
determine the effective behavior of the LCC. 

The emphasis in this paper, however, is on the results rather than the methodology. For this reason, 
the details of the 'second-order' method, which has already been discussed in some detail in the context of 
composites with random microstructures (Ponte Castaheda, 2002a; Idiart and Ponte Castaheda, 2007a), are 
summarized together with the resulting estimates in Appendix A. In view of the high accuracy of the linear 
estimates utilized, any discrepancies between the 'second-order' estimates and the numerical simulations 
observed in Sections 3 and 4 below, in the range < / < 0.3, can be mostly attributed to the error 
introduced in the linearization. 

Finally, it is noted that the strain and stress formulations of the 'second-order' method -based on 
linearizations of (1) and (2), respectively - deliver different predictions: the estimates exhibit a duality gap 
(see Ponte Castaheda (2002a); Idiart et al. (2006a)). Thus, two different sets of 'second-order' estimates will 
be reported, based on the strain (W) and stress (U) formulations. 

2.5. Numerical Fast Fourier Transform method 

The FFT approach used in this work is well-known (Michel et al., 2001; Moulinec and Suquet, 2003, 
2004). For this reason, we shall limit ourselves here to a few remarks. The need to deal with cases of infinite 
contrast and strongly nonlinear behavior (small or zero nonlinearity exponent, m ^ 1) prompted us, mainly 
for convenience, to consider the modified version of the augmented Lagrangian method (Michel et al., 2001) 
proposed by Willot and Pellegrini (2008). 

The L X L pixel medium is discretized with resolutions L = 256, 512 and 1024. The incompressibility of 
the material is approached by a bulk compressibility modulus k = 1000, with elastic shear modulus /i = 1. 
The plastic yield stress was <Jo = 1. Loading is carried out by increasing the driving overall deformation in 
steps of approximately A(s||) = 0.1, up to (sy) = 2 to 4. Such high values were found necessary to obtain 
reliable averages and second moments of the strain. The quantity 77^ = (||V • cr|p)/a : a ^ is used to 
monitor convergence at each loading step, where iterations are stopped for r] of order 10~^ (or even 10~^ 
for the smallest inclusion volume fractions in the stiffest cases) and/or when no notable progress occurs in 
the decrease of 77. In the process, records are made of the iteration number n, and of the fields of interest, 
x^'^\ From three such records (including the most converged one), theoretical converged limits x^^^ are 
estimated by a power-law extrapolation of the type x'^^^L) = x^'^^L) -h a{L)n~^^^^ from which a, /5, and 
x^^^ are determined. A similar extrapolation procedure on L is then used to deduce the infinite-resolution 
limit x^^^ (L ^ 00) from the converged results at the three resolutions. 

The modified version of the method alluded to above, detailed in Willot and Pellegrini (2008), consists 
in replacing the continuum Green function by a discrete counterpart. The discretization scheme enforces 
compatibility and equilibrium at each pixel, consistently with the system discreteness. Probably for this 
reason, the modified algorithm was observed to enjoy faster convergence rates, especially for the stiffest 
cases^. No specific investigation of convergence properties was carried out in the present work, but the 
typical number of iterations needed here per loading step is of order 10^, except when m <C 1 where it can 
climb up to 10^. It is noted however that in terms of CPU, solving problems with m > is much more 
time-consuming than with m = 0. We refer the reader to Eq. (42) in Michel et al. (2001), which has no 
analytical solution for m > 0. As a consequence, it has to be solved numerically at each material point and 
each iteration, whereas an analytical solution is directly implemented for perfectly plastic materials. 

A few cases were checked by comparing to calculations using the continuous Green function. Quite gen- 
erally, excellent agreement was found between both methods save in one remarkable case: in the perfectly- 
plastic limit with rigid inclusions, markedly different solutions are obtained. This difference in behavior 
presumably originates from the short-range spatially-dispersive ("smoothing") features of the discretized 



^For instance, with the above convergence criterion, and in the more demanding situation of a pixelwise disordered medium, 
as was considered in Willot and Pellegrini (2008), the number of iterations to convergence is roughly proportional to 77"°-'^ 
with the discretized Green function, and to ri~^ with the continuous one. 
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Green function (Willot and Pellegrini, 2008), though the precise mechanism at play remains to been inves- 
tigated. More importantly, the solutions from the modified approach turn out to be the most relevant ones 
to the present theoretical study. This point is further discussed in Sec. 4.1. 

3. Results for porous materials 

3.1. Field maps 

Full-field distributions in porous materials under pure shear (PS) and simple shear (SS) loadings have 
been obtained by the FFT method described in Sec. 2.5. For porosities / < 0.3, the simulations show that 
the reduced displacement u* = u — s • x within the unit cell organizes in a periodic array of counter-rotating 
vortices or convection cells. For the case of a linear matrix (m = 1), the patterns can be found in Table I of 
Willot et al. (2008a) (maps B and E). As the nonlinearity increases, the patterns remain similar in character 
but exhibit stronger displacement gradients, and therefore more localized strains. 

The corresponding strain and stress distributions are shown in Table 2, respectively normalized by the 
equivalent macroscopic strain and stress ag. Maps of the 'parallel' (||) and 'perpendicular' (±) components 
(see end of Section 2.2 for definition) are displayed for / = 0.1 and several values of the strain-rate sensitivity. 
Each map is labelled by a number (1-20), and goes along its own field scale on its right, in correspondance with 
the color scale at the extreme right of the rows. As anticipated above, the distributions of both components 
of the strain (maps 1-10) exhibit smooth variations in the linear case but become progressively more localized 
in bands as the nonlinearity increases. Across such bands, the tangential component of the displacement 
field varies significantly. More precisely, if the vectors t and n denote, respectively, the directions tangential 
and normal to the band, then Stn increases with decreasing band width, but not {Snn — ^tt)/2. Indeed, the 
maps show that the deformation tends to localize in bands running across the specimen along directions 
that are parallel to those of maximum macroscopic shear: 0° — 90° and ±45° with respect to the reference 
axes of Fig. 1, under SS (maps 2 and 7) and PS (maps 4 and 9) loading, respectively. As will be seen in the 
next subsection, the fact that these bands have some preferential orientation dictated by the macroscopic 
loading leads to a strong anisotropy of the strain fiuctuations with increasing nonlinearity. 

In the strongly nonlinear ideally plastic limit (m ^ 0), the variation of the tangential displacement across 
the localization bands becomes discontinuous, and consequently the 'parallel' strain becomes unbounded 
along well-defined straight shear bands (maps 1 and 5). The paths followed by these shear bands are such 
that the macroscopic (dissipation) strain energy is minimal (Duxbury et al., 2006). However, equilibrium in 
ideally plastic materials requires that shear bands meet free boundaries at ±45° with respect to the normal 
to the boundary (Kachanov, 2004). Consequently, the bands have to bend in regions sufficiently close to 
the boundary of the pores (see insert in map 1), which causes the 'perpendicular' component of the strain 
(maps 6 and 10) to become very large in those regions^. A closer look at the shear bands reveals a certain 
structure of finite width set by the diameter of the pores, but further work is required to estimate this 
width. Interestingly, while in the SS case the two orthogonal families of bands intersect each other 'inside' 
the pores (map 1), in the PS case they do so at the corners of the unit cell (map 5). These are points 
where the 'parallel' strain takes the highest values. In fact, these are points in the matrix phase where four 
displacement vortices meet, leading to large and highly localized displacement gradients, which are strongly 
enhanced by the constitutive nonlinearity. 

On the other hand, the stress field (maps 11-20) exhibits smooth variations for all m. In the ideally 
plastic limit, the 'parallel' component of the stress (maps 11, 15) takes the highest values in the regions 
where the above-mentioned shear bands develop, as expected. 

Finally, it is noted that the 'minimal paths' followed by shear bands in PS change their configuration for 
porosities / > fc ~ 0.52, see maps 1-6 in Table 3. While for f < fc the two families of shear bands intersect 
at the corners of the unit cell, for f > fc they intersect at the mid-edge points of the unit cell. This is 
because the (dissipation) strain energy, and consequently the total length of the shear bands, must go to zero 



^In the simulations, the strain values within the shear bands are large but finite; however, these values are found to increase 
with increasing mesh resolution, which indicates that the exact solution involves unbounded strains. 
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m = ni = 0.2 rn — 1 m = 0.2 rn = 

Porous - Simple shear | Porous - Pure shear 




Rigid - Simple shear | Rigid - Pure shear 




Table 2: Maps of the 'parallel' (||) and 'perpendicular' (±) components of the shear strain and stress fields, in porous and 
rigidly-reinforced power-law materials under simple shear and pure shear loadings, for / = 0.1 and several values of the strain- 
rate sensitivity (m = 0, 0.2, 1). The strain and stress fields are normalized by the equivalent macroscopic strain £e and stress 
ae, respectively. 

as / approaches the close-packing threshold fcp. Only the latter configuration satisfies this requirement. As 
will be seen in the next subsection, this change in shear band pattern is manifested at the macroscopic level 
as an inflexion point (or kink) in the effective flow stress for PS as a function of / (see Fig. 3a). 
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Table 3: Change of shear band patterns with concentration / in porous and rigidly-reinforced ideally plastic materials under 
pure shear (PS) loading. The strain field is normalized by the equivalent macroscopic strain Se- 



3.2. Effective behavior and field statistics 

The effect of strain-rate sensitivity m on the 'second-order' (SO) and FFT predictions is explored in 
Fig. 2, for a moderate porosity (/ = 0.1). Parts (a) and (b) show plots for the effective flow stress ao 
under 'pure shear' (PS) and 'simple shear' (SS) loadings, respectively, normalized by the flow stress of the 
matrix ctq. As anticipated in Section 2.4, the strain (W) and stress (U) versions of the SO estimates, as 
given by expressions (46) and (72), are not identical; however, this duality gap is very small for most values 
of m. Good agreement is observed between the SO estimates and the numerical simulations for all values 
of m. The predictions are softer for SS than under PS loading. This anisotropy increases with increasing 
nonlinearity and is entirely due to the anisotropic arrangement of the phases: it reflects the fact that most 
of the energy is stored (dissipated) along the localization bands described above, which are shorter in SS 
(see maps 1 and 5 in Table 2). The SO estimates are able to capture this interplay between geometrical 
arrangement and constitutive nonlinearity of the phases in the effective behavior, in contrast to the classical 
Voigt bound, which depends only on the volume fraction of the phases. 

A key improvement in the 'second-order' method over earlier 'linear comparison' methods is the use 
in the linearization scheme of generalized secant moduli whose anisotropy can depend on microstructural 
features and loading conditions, in addition to constitutive nonlinearity. Fig. 2c shows the anisotropy ratios 
k associated with the SO estimates for SS, which solve equations (45) and (71). (The corresponding plots 
for PS are similar and therefore omitted.) Indeed, a non-trivial dependence of /c on m is observed, such 
that the anisotropy strength is increasingly stronger with increasing nonlinearity, the compliant mode being 
set by the direction of macroscopic load {k < 1). In contrast, linearization schemes based on the secant 
and tangent moduli restrict the anisotropy to be /c = 1 and /c = m, respectively, and as a result lead to 
predictions that are less accurate in general, especially for field statistics (see comparisons in Idiart et al. 
(2006b) for random composites). 
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(d) (e) (f) 

Figure 2: Comparisons between the 'second-order' (SO) estimates and full-field FFT numerical calculations for porous, power- 
law materials vs. strain-rate sensitivity m: effective fiow stress ao, normalized by the fiow stress of the matrix phase ctq, for (a) 
'pure shear' (PS) and (b) 'simple shear' (SS) loadings; (c) anisotropy ratio k of linear comparison matrix associated with SO 
estimates for SS; (d)-(f) field statistics for SS, normalized by the macroscopic equivalent strain Se and stress ae- The porosity 
is / = 0.1. 

Predictions for the field statistics under SS loading are shown in Figs. 2d-f. (The corresponding results for 
PS are similar and therefore omitted.) In these plots, the strain and stress quantities are normalized by the 
macroscopic equivalent strain Se and stress ag, respectively. The FFT results show that as the nonlinearity 
increases, the pores undergo larger deformations and the strain fluctuations in the matrix increase, while the 
stress fluctuations remain fairly constant. The strong anisotropic growth of strain fluctuations is a direct 
consequence of the strain localizing along bands with a preferred orientation determined by the loading, as 
observed in Table 2 (maps 1-10). In the ideally plastic limit, the strain fluctuations become unbounded"^. 
The SO estimates, given by expressions (38)- (39) and (67)- (68), are in good agreement with the numerical 
results in the range 0.3 ^ m < 1. Within that range, the anisotropy ratios of the LCC matrix lie in the range 
0.3 ^ /c < 1, see Fig. 2c, and the local fields in the LCCs, reported in Willot et al. (2008a)^, mimic fairly 
well the nonlinear fields. For smaller values of m, however, the linear fields cannot mimic as accurately the 
localized character of the nonlinear strain field, and the SO predictions for the field fluctuations deteriorate 
significantly. In particular, the SO predictions for the strain fluctuations remain bounded as m ^ 0. We 
note in passing that strain statistics are more accurately predicted by the strain-based SO(W) estimates, 
while stress statistics are more accurately predicted by the stress-based SO(U) estimates. 



^The FFT results for these quantities are finite, but increase with increasing mesh resolution. 

^ Table II in Willot et al. (2008a) shows stress maps for various values of k; due to linearity, strain maps are the same up to 
a change of scale given by the matrix moduli. 

10 



0.8 



0.6 



0.4 



0.2 







% Voigt 








\ V ■• LA-UB 






SO(U) 




FFT 


PS LA-LB 




m = 





0.2 0.4 0.6 0.8 

/ 



(a) 




Figure 3: Comparisons between the 'second-order' (SO) estimates and full-field FFT numerical calculations for porous, ideally 
plastic materials vs. porosity /: effective flow stress ao for (a) 'pure shear' (PS) and (b) 'simple shear' (SS) loadings; (c) 
anisotropy ratio k of linear comparison matrix associated with SO estimates for PS. 
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Figure 4: Comparisons between the 'second-order' (SO) estimates and full-field FFT numerical simulations for power-law 
materials with dilute inclusion concentration: (a) porous materials and (b) rigidly-reinforced composites. Also included for 
comparison are the results of Fleck and Hutchinson (1986) and Lee and Mear (1992). 



The effect of porosity in ideally plastic solids (m = 0) is explored in Fig. 3. Parts (a) and (b) show plots 
for ao, and include the rigorous upper (UB) and lower (LB) bounds derived in Appendix B by limit analysis 
(LA). The FFT results lie close to the LA upper bound, and are consistent with a close-packing threshold 
given by fcp = ir/A ^ 0.785, the porosity at which neighboring pores come into contact with each other. 
As anticipated in the previous subsection, the FFT results for PS exhibit an inflexion point (or a kink) at 
/c ~ 0.52, as a consequence of a change in the shear band pattern, see Table 3. The LA upper bound given 
by (80) also exhibits this feature. It is quite remarkable that both versions of the SO estimates are able to 
capture, at least qualitatively, the correct anisotropic dependence of 5o on porosity, including the inflexion 
point in the PS case. At large porosities, however, both SO estimates become rather stiff and remain finite 
at the close-packing threshold. This inaccuracy is inherited from the linear estimates used for the LCC, 
which become inadequate at porosities close to fcp (see Willot et al. (2008a)). Improved SO estimates for 
large porosities could thus be obtained by using more appropriate/accurate linear estimates. 

3.3. Dilute limit 

Consider power-law materials with vanishingly small porosity levels (/ 0) and fixed rate-sensitivity 
m > 0. For these materials, it is widely accepted that the effective behavior is insensitive to the actual 
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distribution of porosity. Following Duva and Hutchinson (1984), Fleck and Hutchinson (1986) proposed the 
ansatz _ 

^^l-a(m)/, (13) 

and determined the coefficient a{m) numerically by considering an isolated circular pore in an infinite 
matrix. The present study confirms the ansatz (13): both the FFT results and the SO estimates exhibit a 
linear correction at vanishingly small porosities. The predicted behavior is thus isotropic. In particular, the 
SO estimate for a{m) can be obtained by expanding (46) and (72), see Appendix A for details; the result is 

(1 + v^)' 

for both versions. A comparison with the numerical results for a is shown in Fig. 4a (Fig. 4b will be 
discussed in Section 4). As can be seen in the figure, the numerical and theoretical results coincide for 
m = 1, agreeing with the classical result of Eshelby (1957), and are in good agreement for all other values 
of m. Both the FFT simulations and the SO estimates predict vanishingly small strain fluctuations in the 
matrix phase as / ^ 0, which is consistent with non-interacting pores in the dilute limit. As m ^ 0, 
however, a{m) becomes unbounded, and the range of validity of (13) vanishes. 

Indeed, for ideally plastic solids (m 0) with small but fixed porosity / > 0, the limit analysis bounds 
derived in Appendix B demand that 

^ = l-ao/^/^ (15) 

where the coefficient ao depends on the specific loading condition (PS or SS), see expressions (79) and (80). 
Thus, the dependence of 5o on f in the dilute limit is now non- analytic^ and the behavior is anisotropic. This 
is a consequence of the developing strain localization, see maps 1-10 in Table 2, which introduces long-range 
interactions between pores. As a result, the effect of the actual distribution of porosity on the effective 
behavior persists in the dilute limit, and a crossover between the regular regime (13) and a singular regime 
(15) ensues. One can attempt to estimate the crossover threshold /c(^) between these regimes by equating 
expressions (13) and (15), to find 

/c(m) = aQ/a{m)'^ ~ Aa^m for m <C 1. (16) 

However, preliminary results for sufficiently small / and m -not shown here- indicate that the crossover 
threshold may be much lower. A detailed study of the crossover will be reported elsewhere. 

The FFT results for dilutely voided ideally plastic materials are consistent with (15), with the coefficient 
Qfo being very close to that of the upper bound - a numerical fit to the results yields o^o = 0.826 and 
Qfo = 1.185 in PS and SS, respectively. The upper bound coefficient is given in terms of the ratio between 
the length of the bands and the side length of the unit cell / by o^o = {'^l - The ratio is Z^/Z = 1 
for SS, and h/l = for PS, see Figs. 7a, b. The porous material is thus weaker under SS than under PS 
because the underlying shear bands are shorter in the ffist case. In turn, the exponent 1/2 in (15) arises 
because the energy is dissipated at one-dimensional straight shear bands running through two-dimensional 
pores (Drucker, 1966). As a result, the effective behavior remains anisotropic in the dilute limit. In addition, 
the strain fluctuations in the matrix phase remain unbounded as / ^ 0, which is consistent with a collective 
behavior of the pores in the dilute limit. 

'Second-order' estimates for dilutely voided ideally plastic materials are obtained by expanding (46) and 
(76), resulting in the expression 

^^l-asof/\ (17) 

where the coefficient is independent of the loading, and is given by ago = 1/2 in the strain version and 
aso = 3/2^/^ in the stress version. In addition, the SO estimates predict non- vanishing strain fluctuations 
in the matrix phase as / ^ 0. Thus, remarkably, the SO estimates also predict a crossover between the 
regular regime (13) and the singular regime (17), although the exponent in the singular regime is incorrect. 
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In this connection, it is recalled that SO estimates for 2D random porous materials exhibit exactly the same 
dilute hmit (Ponte Castaheda, 2002b; Pastor and Ponte Castaheda, 2002). The 2/3 exponent thus appears 
to be independent of the distribution of porosity. 

It is noted, however, that the linearization in both versions of the SO estimates is such that A: ^ 
as / ^ 0, see insert in Fig. 3c. In this limit, the LCC effective moduli (40) exhibit a crossover behavior 
(see Willot et al. (2008a)): if ^ tt/, the medium is strongly anisotropic, the strain in this linear medium 
localizes in bands connecting the pores, and the dilute limit depends on the distribution of porosity; if 
on the other hand k > nf, the anisotropy is weak and the pores behave as isolated inclusions. Making 
use of the asymptotic expansions of Willot et al. (2008a), it is shown in Appendix A that k ^ /^/^ for 
both versions. Thus, the dilute SO estimates (17) do not depend on the actual distribution of porosity 
because they make use of the 'uncor related' regime of the linear estimates. However, the fact that k ^ 
suggests that the 'correlated' regime could be exploited by suitably modifying the linearization scheme. In 
this connection, it should be noted that the estimates (17) make use of an ad- hoc prescription for certain 
reference tensors entering in the linearization, see expressions (32) and (61). The above findings hint that 
alternative prescriptions may be developed resulting in dilute SO estimates that depend on the distribution 
of porosity (in the singular regime), and that would yield improved estimates more generally. 

4. Results for rigidly-reinforced composites 
4.1. Field maps 

Full-field distributions in rigidly-reinforced materials are shown in Table 2 (maps 21-40), for three values 
of the strain-rate sensitivity {m = 1,0.2,0) and a moderate reinforcement concentration (/ = 0.1). Unlike 
in the porous case, the strain does not localize in shear bands as the nonlinearity increases. Instead, an 
inclusion screening effect develops with increasing nonlinearity: the strain becomes small (compared to €e) 
along bands of width equal to one inclusion diameter, running across the specimen along the directions 
of maximum macroscopic shear (0° and 90° with respect to the Cartesian axes) and seeking out the rigid 
inclusions, see maps 21-25 in Table 2. Thus, the macroscopic deformation is mainly accommodated by a fairly 
uniform strain present in the remaining regions outside these bands. The corresponding stress distributions, 
on the other hand, are seen to become progressively more uniform as the nonlinearity increases, see maps 
31-40. In the ideally plastic limit, both components of the stress are uniform almost everywhere, exhibiting 
variations only in small regions around the inclusions (maps 31, 35, 36, 40). In fact, the matrix phase is at 
yield almost everywhere, and this translates into no reinforcement effect in the effective behavior due to the 
rigid inclusions, as will be seen below. 

Under PS loading, the two families of bands intersect each other not only at inclusion sites but also 
at regions in the matrix phase (corners of the unit cell), giving rise to the formation of square regions 
where the strain is almost null (see map 25). No such regions develop under SS shear (see map 21). This 
difference is responsible for the strong anisotropy in the effective behavior observed in Fig. 5 below for high 
reinforcement concentrations. Indeed, while in the SS case the strain pattern remains essentially unchanged 
for any reinforcement concentration below the close-packing threshold fcp = 7r/4 ^ 0.78, in the PS case 
the strain pattern has to change abruptly for reinforcement concentrations larger than / = 7r/8?^0.39, the 
value of / at which the square regions come into contact with each other. This transition is shown in Table 
3, maps 7-12. For / > 0.4, the strain tends to concentrate in a rim around the inclusions, being very small 
in most of the matrix phase. As will be seen below, this results in a significant increase in the effective fiow 
stress for that range of reinforcement concentrations. See also Willot et al. (2008b) for analogous effects in 
the linear anisotropic porous case. 

It should be emphasized at this point that, due to the lack of strict convexity of the potentials (1) and 
(2) when m = 0, the strain field in an ideally plastic composite may not be unique (Suquet, 1981) (the 
effective stress, however, is uniquely defined, Bouchitte and Suquet (1991)). In fact, it is easy to verify that, 
since the matrix of the reinforced materials considered here is at yield almost everywhere, alternative strain 
fields to the 'diffuse' fields shown in Table 2 are possible and would consist of infinitely thin straight shear 
bands running through the matrix at 0°-90° for SS, and ±45° for PS (provided / < n/S). In this connection. 



13 



1.3 
1.25 

1.2 
1.15 

1.1 
1.05 















PC 






f=o.i 


SO(U) 






SO(W) - 


FFT 




1 




Reuss 


t 



0.2 0.4 0.6 0.^ 

(a) 



2.5 



1.5 



0.5 



SS 
















SO(W) 




SO(U) 




. FFT 


/=0.1 



0.2 0.4 0.6 0.^ 

(d) 



1.3 
1.25 

1.2 
1.15 

1.1 
1.05 











- cro 






0.8 








0.6 


f=o.i 


SO(U) 


SO(Wl.^- 


0.4 


FFT 


Reuss 


'"o 


0.2 



0.2 0.4 0.6 

(b) 



0.8 1 
m 



0.4 
0.35 

0.3 
0.25 

0.2 
0.15 

0.1 
0.05 




: SS 




SO(W} 






SO(U) 




• 


FFT : 




SD'i'C^II )..... 














f=0.l 



0.2 0.4 0.6 0.^ 

(e) 



0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 




k 




SO(W) // 








// SO(U) 






SS 
/=0.1 


) 0.2 0.4 0.6 

(c) 


0.8 ] 
m 


SS 


SO(W) 

SO(U) 

. FFT 










0.2 0.4 0.6 


0.8 



(f) 



Figure 5: Comparisons between the 'second-order' (SO) estimates and full-field FFT numerical calculations for rigidly- 
reinforced, power-law materials vs. strain-rate sensitivity m: effective fiow stress ?o for (a) 'pure shear' (PS) and (b) 'simple 
shear' (SS) loadings; (c) anisotropy ratio k of linear comparison matrix associated with SO estimates for SS; (d)-(f) field 
statistics for SS. The reinforcement concentration is / = 0.1. 



it was found that the same augmented Lagrangian method, separately apphed with the discretized Green 
function (see Section 2.5), and with the 'continuous' Green function (as used in the original FFT algorithm 
of Michel et al. (2001)) converge to completely different 'diffuse' and localized solutions for the strain field, 
respectively (whereas the stress field, is independent of the choice of the Green function at high resolution). 
With the 'continuous' green function, thin straight shear bands tangent to the inclusion-matrix interface 
are observed (see also Moulinec and Suquet (1998)). In any event, as will be seen below, the ideally plastic 
'second-order' estimates derived in this work, which follow from a power-law regularization, predict strain 
statistics entirely consistent with the 'diffuse' strain fields of Table 2, and not with the localized solution. 

4.2. Effective behavior and field statistics 

The effect of matrix nonlinearity m on the 'second-order' (SO) and FFT predictions is explored in Fig. 
5 for a moderate reinforcement concentration (/ = 0.1). Parts (a) and (b) show plots for the effective flow 
stress ao under 'pure shear' (PS) and 'simple shear' (SS) loadings, respectively, normalized by the flow stress 
of the matrix ctq. Both versions of the SO estimates are found to be in good agreement with the FFT results, 
for all values of the strain-rate sensitivity m and both loading conditions. The estimates agree with the FFT 
results in that ao decreases monotonically with increasing nonlinearity (i.e., decreasing m), and in that the 
composite material is 'stronger' under PS than under SS loadings, for all values of m different than zero. 

In the ideally plastic limit (m 0), the limit analysis upper bound derived in Appendix B coincides 
with the classical Reuss lower bound, and therefore give the exact result which is ao/c^o = 1 for both loading 
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conditions. Indeed, this is the result obtained with the FFT calculations. Both versions of the SO estimates 
are found to coincide with the exact result in this strongly nonlinear limit, thus predicting no reinforcement 
effect due the presence of rigid inclusions. In this connection, it is recalled that this is the exact result 
whenever the macroscopic deformation can be accommodated by straight shear bands passing through the 
matrix phase and avoiding the inclusions (see Drucker (1966) and Taliercio (1992)). In the SS case, this is 
possible for any reinforcement concentration smaller than the close-packing threshold /cp, while in the PS 
case this is not possible for / > tt/S and therefore a finite reinforcement effect is expected in that case. 

Predictions for the field statistics under SS loading are shown in Figs. 5e-f, normalized by Se and ag. 
Unlike in the porous case, the SO estimates are in good agreement with the FFT results for all values of 
m. This is mainly because the strain field in the LCCs can now mimic fairly well the nonlinear strain field, 
even in the ideally plastic limit. Indeed, as m ^ the anisotropy ratio /c ^ 0, see Fig. 5c, and the strain 
distributions in such LCCs, reported in Willot et al. (2008a)^, are similar to those of maps 21-30 in Table 2. 
The nonlinear predictions show decreasing inclusion stress concentration and matrix stress fluctuations with 
increasing nonlinearity. This is consistent with the stress field becoming more homogeneous throughout the 
composite as m ^ 0, see maps 31-40 in Table 2. The small discrepancies between the SO estimates and the 
FFT results in the ideally plastic limit are due to the small local variations of stress around the inclusions 
observed in the numerical simulations. The predictions also show that the strain fluctuations remain bounded 
in the ideally plastic limit, unlike in porous materials. This is because the strain distribution is essentially 
piecewise uniform rather than localized in shear bands (cf. maps 21-30 and maps 1-10 in Table 2). In this 
connection, it is interesting to note that recent FFT simulations of rigidly-reinforced composites with random 
microstructures (Moulinec and Suquet, 2003; Idiart et al., 2006b) have shown shear band strain localization 
and unbounded 'parallel' strain fluctuations as m ^ 0, in agreement with the 'second-order' predictions for 
that class of materials (Idiart et al., 2006b; Idiart and Ponte Castaheda, 2007b). It is thus remarkable that 
the 'second-order' method, being only a homogenization theory, is able to correctly predict the completely 
different character of strain distribution arising in periodic versus random reinforced composites, as reflected 
not in the macroscopic behavior but rather in higher-order statistical information. 

The effect of reinforcement concentration is explored in Fig. 6, for a moderate nonlinearity (m = 0.2). 
Parts (a) and (b) show plots for ao- The FFT results show that the composite becomes 'stiffer' as the 
reinforcement concentration increases, as expected, and that, for moderate and higher concentrations, this 
effect becomes significantly stronger for PS than for SS loadings. As anticipated in the previous subsection, 
this nonlinearity-induced anisotropy in the effective behavior is due to the fact that, as the nonlinearity 
increases, the circular inclusions behave progressively more like square inclusions with an orientation relative 
to the lattice that depends on the loading orientation (see Table 2). Due to the orientation of these 'square' 
inclusions, the close-packing effects become important at lower values of / in PS than in SS loadings. 
This, together with the fact that the linear estimates utilized in the context of the SO estimates become 
progressively less accurate as the close-packing effects become more important (see Willot et al. (2008a)), 
is the reason why, at large values of /, the SO estimates are significantly less accurate for PS than for SS 
loadings. As already mentioned in the previous section, the SO estimates could be improved for high values 
of / by using a different linear estimate, appropriate for reinforcement concentrations near the close-packing 
threshold. 



J1..3. Dilute limit 

Consider power-law materials with vanishingly small reinforcement concentrations (/ ^ 0) and fixed 
strain-rate sensitivity m > 0. As already discussed, for these materials, the effective behavior is expected to 
be insensitive to the actual distribution of inclusions, and Lee and Mear (1992) proposed the ansatz 

^«l + /3(m)/, (18) 



^ Table II of Willot et al. (2008a) shows fields for porous materials with a different definition for fc, and require suitable 
reinterpretation. Thus, maps 3-5 and 8-10 in that table can be interpreted as strain fields in rigidly-reinforced composites 
under PS, while maps 16-18 and 21-23 can be interpreted as corresponding strain fields under SS. 
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Figure 6: Comparisons between the 'second-order' (SO) estimates and full-field FFT numerical calculations for rigidly- 
reinforced, power-law materials vs. reinforcement concentration /: effective flow stress ?o for (a) 'pure shear' (PS) and 
(b) 'simple shear' (SS) loadings; (c) anisotropy ratio k of linear comparison matrix associated with SO estimates for PS. The 
strain-rate sensitivity is m = 0.2. 



and determined the coefficient P{m) numerically by considering an isolated circular inclusion in an infinite 
matrix. The present study confirms the ansatz (18): both the FFT results and the SO estimates exhibit 
a linear correction at vanishingly small /. The predicted behavior is thus isotropic. In particular, the SO 
estimate for /3{m) can be obtained by expanding (46) and (72), noting that A: ^ m as / ^ (see Fig. 6c); 
the result is 

P{m) = iy^(l + V^)2 (19) 

for both versions. This result is in exact agreement with that obtained in the context of random rigidly- 
reinforced composites (Ponte Castaheda, 1996, 2002b). A comparison between the SO estimate and the 
numerical results for f3 is given in Fig. 4b. The numerical and theoretical results coincide for m = 1 with 
the classical result of Eshelby (1957), as they should, and are in excellent agreement for all other values of 
m. Note that, unlike in the porous case, the coefficient (19) vanishes as m ^ 0, suggesting an extended 
range of validity of (18) in this limit. 

Now consider ideally plastic solids (m 0) with fixed reinforcement concentration 0</<7r/8. As 
already mentioned in Section 4.2, the limit analysis bounds derived in Appendix B demand that 

^ = 1, (20) 

independent of the value of /, in exact agreement with the corresponding FFT results and SO estimates. The 
dilute expansion of (20) is consistent with the ideally plastic limit of (18). Thus, in contrast to the porous 
case, the dilute concentration and small rate-sensitivity regimes in reinforced composites are consistent, at 
least to first order in /. In addition, the FFT simulations and SO estimates predict vanishingly small strain 
fiuct nations in the matrix phase as / ^ for both regimes. This suggests that, in dilutely reinforced ideally 
plastic solids, only weak long-range interactions between inclusions develop, and the inclusions effectively 
behave as if they were isolated at leading order. This in turn may explain why the 'second-order' predictions 
are more accurate for reinforced composites than for porous materials. It is noted, however, that the 
present analysis does not exclude the existence of a crossover behavior at higher-order corrections. In fact, 
preliminary results for sufficiently small / and m do suggest that such a crossover exists, and will be reported 
elsewhere. 



5. Concluding remarks 

This paper presented a combined numerical-theoretical study of the macroscopic behavior and local 
field distributions for a special class of two-dimensional periodic composites with viscoplastic phases, with 
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emphasis on the infinite contrast, strongly nonlinear case. The comparisons show good overall agreement 
between the 'second-order' (SO) homogenization estimates and the FFT full-field numerical simulations, not 
only for the effective behavior but also for the first- and second-order field statistics of the stress and strain 
fields. 

More specifically, the strain field for the rigidly-reinforced materials was found to develop (with decreasing 
strain-rate sensitivity m) a banding pattern consisting of low-strained bands of finite width set by the size 
of the inclusions. For both pure and simple shear loading, and volume fractions allowing the passage 
of straight bands through the matrix, the resulting interaction between inclusions is not strong, and the 
effective fiow stress tends to the yield stress of the matrix phase in the ideally plastic limit. In these cases, 
the SO estimates are extremely accurate. However, for pure shear and volume fractions greater than 7r/8, 
the inclusions block the straight paths and the deformation mechanism changes abruptly leading to a finite 
strengthening effect for larger volume fractions all the way up to inclusion contact (which the SO estimates 
only capture approximately). By contrast, the strain field in porous materials tends to localize (again with 
decreasing strain-rate sensitivity) in bands seeking to pass through the pores in order to minimize the plastic 
dissipation. In this case, the interaction between the pores is obviously quite strong and dependent on the 
geometry of the configuration as well as the relative orientation of the loading (e.g., pure versus simple 
shear). Again for pure shear, an abrupt change in deformation mechanism is observed for values of the 
porosity near 50%. In this case, the SO estimates are quite accurate for values of strain-rate sensitivity 
between 0.2 and 1, and deteriorate for smaller strain-rate sensitivities (although they correctly capture the 
transition). 

One of the main findings of this work is the observation of a crossover-type behavior for the porous case 
in the dilute concentration (/) and small rate-sensitivity (m) regimes. Indeed, the FFT simulations show 
that for finite values of m > 0, the effective fiow stress relative to the flow stress of the matrix is reduced from 
unity by a factor proportional to the porosity / (cf. (13)). This is consistent with the classical dilute limit 
with non-interacting inhomogeneities. However, the coefficient a (see (14)) blows up as m~-^/^ when m ^ 0, 
suggesting that the range of validity of this expansion tends to zero in the the ideally plastic limit. In fact, 
first taking the limit as m ^ 0, and then considering small porosity / leads to a non- analytic dependence on 
the porosity (1/2-power dependence on f - see (15)), together with non- vanishing strain fluctuations in the 
matrix phase, which is consistent with a collective-type behavior induced by the strong interactions among 
the voids as the developing shear bands seek to minimize the overall plastic dissipation by passing through 
the voids. Between these two asymptotic regimes, there is an 'intermediate' asymptotic (i.e., crossover) 
regime where the type of solution transitions from one regime (i.e., 'strongly dilute') to the other (i.e., 
'strongly nonlinear'). Mathematically, this phenomenon can be linked to the change of character of the 
governing equations — from elliptic to hyperbolic (see Chenchiah and Bhattacharya (2005) and Willot et al. 
(2008b) for analogous results in other systems). In this connection, it should be mentioned that, remarkably, 
the SO theory is able to predict a crossover behavior for the porous materials, with a very good quantitative 
prediction in the strongly dilute regime (see Fig. 4a for ce), but with the wrong exponent (2/3 versus 1/2) 
in the strongly nonlinear regime. However, the true form of the crossover in this case is still unknown. 

For the rigidly-reinforced materials, on the other hand, the two asymptotic regimes were found to be 
consistent up to leading order. In this case, the strongly dilute regime is linear in / (cf. (18)) with a 
coefficient P{m) that tends to zero as m ^ 0, while the strongly nonlinear regime has a vanishing correction 
for sufficiently small values of / (cf. (20)). In addition, the strain fiuctuations are vanishingly small in this 
limit. This suggests that the ideally plastic limit in the case of rigid inclusions is in some sense less singular, 
which may help explain why the SO theory is able to deliver much more accurate predictions in this case. 
We note, however, that the analysis provided in this work does not exclude the possibility of a crossover- type 
behavior in the higher-order corrections, which will be investigated in future work. 

It should be emphasized that the combined dilute and ideally plastic limit for porous media is particularly 
relevant to the ductile fracture initiation problem, where the distribution of porosity is expected to be 
random. However, for random microstructures, the corresponding numerical simulations are much more 
challenging than for the periodic microstructures discussed in this work. The only results that the authors 
are aware of for porous materials with random microstructures are the numerical results (based on limit 
analysis) of Pastor and Ponte Castaheda (2002), which suggest that the dilute exponent in the 'strongly 
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nonlinear' regime should be between 1/2 and 2/3, at least for a class of 2D microstructures known as 
composite cylinders, as well as the FFT results of Willot (2007) for 'pixel-type' 2D random microstructures 
that suggest an exponent of 2/3. On the other hand, the SO estimate of the Hashin-Shtrikman type predicts 
(Ponte Castaheda, 2002b) a crossover behavior with an exponent of 2/3 in the strongly nonlinear regime 
for 2D random microstructures. Although the role of disorder on the dilute exponent still remains to be 
investigated in detail, the relative success — at least in qualitative terms — of the SO theory for the periodic 
case suggests that the corresponding predictions of Ponte Castaheda (2002b) for the random case should at 
least be qualitatively correct. 

Finally, it should be noted that this study was restricted to isochoric loadings 'aligned' with the periodic 
microstructure. It is emphasized, however, that the homogenization methods utilized can handle completely 
general loading conditions (see Bilger et al. (2005); Danas et al. (2008b)). In particular, porous materials 
subjected to non-isochoric loadings, where more complicated localization patterns arise, will be considered 
elsewhere. 
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A. Second-order estimates for power-law porous materials 
A.l. Strain formulation 

When specialized to porous materials with an isotropic matrix phase characterized by (1), the 'second- 
order' variational estimate for the effective strain potential is given by (Ponte Castaheda, 2002a; Idiart et 
al., 2006a) 

w(e) = Stat \wL (e; e^^\ L^^) ) + (1 - f)v^^^ {e^^\ L^^) ) | , (21) 

where the stationary operation consists in setting the partial derivative of the argument with respect to 
the variable equal to zero. In this expression, -2;^ is the effective potential of a porous LCC with the same 
microstructure as the nonlinear material and with a matrix characterized by a second-order, Taylor- type 
expansion of the nonlinear potential (1)2 about a reference strain tensor e^^\ given by 

£W,L«) = «;«(£W) + ^{e(^)):(e - e^'^) + ^(e - £W):LW :(£ - (22) 

OS z 

where L^^^ is an incompressible, symmetric, fourth-order tensor (of moduli) of the form 

L^^) = 2AE|| +2/iE_L, k = X/fi. (23) 
Here, A and /i are two shear moduli, and 

2 s^^^ s^^^ 

^ll = 3Tfe^Tfe' IE±=IK-E||, (24) 

are two orthogonal, fourth-order, projection tensors with principal axes 'aligned' with £^^\ and K denotes 
the standard, fourth-order, shear projection tensor. It should be noted that, even though the nonlinear 
potential (l)i is isotropic, the tensor of moduli (23) is generally anisotropic. The ratio k constitutes a 
measure of this anisotropy. 

In turn, the 'error' function v^^"^ in (21), a measure of the degree of nonlinearity of w^^\ is defined as 



vW (£-(1) , L(i) ) =stat {^(^) (f^^)) -4'^ (f^^) ; e^'\ L^^)) } 
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(25) 



The nonlinear estimate (21) requires a linear estimate for the effective potential wl of the porous LCC. 
Noting that the matrix potential (22) corresponds to a 'thermoelastic' material with a 'thermal stress' and 
a 'specific heat' given by 



(26) 



T« = ^(e«)-L«:e«, 
g(i) = «;(i)(£(i))-TW:£W-i£W:LW:£W (27) 



the effective potential becomes (Willis, 1981) 

WL{e) = ^6:L:6 + r:6 + 5^, (28) 

with r and g given in terms of L by 

r = L:L(i)-i:r(^\ (29) 
g = (l-/)^a)-r^i):L(i)-i:[L-(l-/)L(i)]:L(i)-i:r(i). (30) 

In these expressions, L denotes the effective modulus tensor of a porous material with a purely elastic matrix 
with modulus tensor (23), and completely characterizes (28). 

For later use, it is recalled that the first and second moments of the strain field Sl in the matrix of the 
porous LCC, can be extracted from the potential wl^Y means of the identities (Idiart and Ponte Castaheda, 
2007a) 



.(1) 



dwr 



1 - f ar(i) 



/ K> \(i) 2 



(31) 

r(i),^(i) 



where the variables that are held fixed have been emphasized. Thus, the variational procedure employed 
to derive the estimate (21) dictates an 'optimal' choice for the moduli A and /i of the matrix phase in the 
LCC, which follows from the stationarity condition in (21). In addition, the reference strain s^^^ must also 
be specified in order to characterize the potential w^j^"^ . Unfortunately, it has not yet been possible to find 
a similar 'optimality' condition for this variable, and an ad hoc prescription for e^^^ is therefore required. 
Hereafter, following Idiart et al. (2006a) we identify 6^^^ with the macroscopic deviatoric strain, i.e., 

£(')=£d. (32) 

In addition to giving sensible results in the case of isotropic composites, this prescription has the advantage 
of simplicity. It is emphasized, however, that this prescription is probably not 'optimal'. 

Then, the stationarity conditions in (21) and (25) yield a system of nonlinear algebraic equations for the 
variables £^^\ A and /i. More explicitly, making use of (1)2, (23) and (32), the conditions resulting from the 
stationary operation in (25) are 

3A (4;) - £e) = (4;V4'^) 4>'{i^^^) - <^'(^e), 3/x= (1/4^)) <^'(4i )), (33) 

where e\P, = \/ (2/3)s*^^^ :E|u :s*^^^ represent the components of s*^^^ that are 'parallel' and 'perpendicular' 



to the macroscopic strain e. The equivalent part of s^^^ is then (fe^^)^ = (^y^^)^ + (f^i^^)^. 
In turn, the conditions resulting from the stationarity operation in (21) are 



.(1) _ / 1 2a^L .(1) / 1 2dwL 



where a choice of roots has been made (Idiart et al., 2006a). The arguments inside the square roots in 
these relations depend on certain traces of the strain covariance tensor in the matrix phase of the LCC 
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(Ponte Castaneda, 2002a) as identities (31) show. It then foUows from relations (33) that the moduh A 
and /i depend on the intraphase field fluctuations (in the LCC). This is one of the main reasons for the 
improvement of the 'second-order' method over earlier linear comparison methods which made use of only 
the first or the second moments of local fields (Ponte Castaneda and Suquet, 1998). 
Finally, making use of the relations (33), the estimate (21) reduces to 



w{e) = (1 - /) [^(ef )) - ^'(ee) {e^^ - e^O] ' (^5) 



where e^^l is the equivalent part of (31) i. 

Effective behavior. An estimate for the effective behavior of the nonlinear porous material is obtained by 
differentiating expression (21) with respect to e. For the specific choice of reference tensor (32), this yields 
(Idiart and Ponte Castaneda, 2007a) 

a=^{e)=WL + {l-f)p'-'\ (36) 
where Wl = Ls + r is the macroscopic stress in the LCC, and the (incompressible) tensor p^^^ is given by 
L« -L(i)1 :(f(i) -eW) + l:^[((£^^-£,)0(£^^-£,))«-(£W-£,)®(fW-£,)l :s,, (37) 



where L^^^ = d'^w^^\e) / dede is the tangent modulus tensor of the matrix phase evaluated at e. The 
averaged quantity appearing in the second term is determined from (31). 

Field statistics. In addition to the effective behavior, estimates for the field statistics can also be extracted 
from corresponding estimates for w by means of the procedure proposed by Idiart and Ponte Castaneda 
(2007a). Applying this procedure in the context of the 'second-order' method, the following estimates for 
the first and second moments of the fields in each phase are obtained (Idiart and Ponte Castaneda, 2007a): 

s(^) =4^\ a^'^ p^'\ (38) 

(S S)(^) = {SL {cT = {CTL CT l)^'^ + 2&^'^ 0, p^^^ , (39) 

where a^^^ = dw^^\e)/de, in the porous phase a^^^ = p^^^ = 0, and the subscript s in the tensor product 
denotes symmetrization. The strain statistics in the LCC are given by expressions (31), and the corre- 
sponding stress statistics follow from the linear stress-strain relation associated with (22). Thus, while the 
nonlinear estimates for the strain statistics in the matrix coincide with those quantities in the LCC, the 
estimates for the stress statistics exhibit certain correction terms which depend on the matrix nonlinearity 
through the tensor (37). 

Expressions for power-law materials under 'aligned' loadings. When the macroscopic loading is 'aligned' in 
the sense described in Sec. 2.1 (see Fig. 1), the (compressible) effective tensor L exhibits the same symmetry 
as the L^i) defined by (23) with (24) and (32), i.e., 

L = 2SJJ + 2AE|| + 2/IE_L, (40) 

where ^ is a bulk modulus, A and Jl are the 'parallel' and 'perpendicular' shear moduli, and JJ denotes the 
standard fourth-order hydrostatic identity tensor. In this work, the moduli are determined approximately 
using the hnear estimates of the Hashin-Shtrikman type derived in Willot et al. (2008a). For SS loadings, 
E||^^ coincide with the tensors E^^'^^ defined by (9), so that the moduli are given by expressions (22) of 
Willot et al. (2008a), with m = ^ = 0; for PS loadings, coincide with the tensors E^^'^^, so that the 

moduli are given by those same expressions with A (A) and /i (Jl) interchanged. 

Evaluating the derivatives of wl in expressions (34) and (31), we arrive at the following relations: 



efV^e = 1 - J (1 + rf f|| + (l - , e^/^e = |1 + r| f'f, = 1 + (1 + r) f, (41) 
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where r = (r^^) ■ T(i))i/V(3Aee), and where 



(1-/)^ 



(1-/) 



A 




(1-/) \A, 



1 A 



(1-/)A 



1. (42) 



Here, the prime denotes differentiation with respect to k. Making use of (33), one obtains r = {1/k 
/se), which, after some algebraic manipulations, yields 



e«/ee = l-r/(l-fc), e^l^/ee = [{l-r)/k]f'f, = 1 + [(1 - r)/fc] f, 

where the function r has been defined as 



r = 



1 + 



k 



1-k 



(43) 



(44) 



Thus, the right hand sides of expressions (43) are explicit in k. Then, by taking the ratio of the two 
conditions (33) with given by (5)i, we obtain the relation 



In 



(1-A:) j^V^e) + A:' 



In 



(45) 



where = (^[i"*"^)^ + i^i^)'^- This relation, together with (43), constitutes a single nonlinear, algebraic 
equation for the anisotropy ratio k that must be solved numerically. (Alternatively, since the right hand 
side of (45) depends on k but not on m, we can assign values to k and obtain from (45) the corresponding 
values of m.) 

Finally, the 'second-order' estimate (35) can be written as (6)i with the effective flow stress given by 



(1-/) 



l+m 



.(1) 



(l + m) 4- 



.(1) 



(46) 



where fj^^ and 6:^^^ are given in terms of k by relations (43). The corresponding estimates for the 

field statistics follow from relations (38)- (39). Since the phase averages should be of the form (11), e'^^^ is 
completely determined by (41), and s^^^ can be readily obtained from the macroscopic balance (1 — f)e^^^ + 
= s. In turn, the estimates for the standard deviations of the 'parallel' and 'perpendicular' components 
of the strain in the matrix phase, as defined by (12), read: 



5D«(£||,x) = |l + r| f\'X. 
In addition, the tensor p^^^ required to compute the stress statistics simplifies to 

p(i)=2(A-A,) (4^)-ei\))(£,/ee), 



(47) 



(48) 



where 2At = Ey • L^^^ = (2/3)(/)'(^e)- Since phase 2 is vacuous, the phase averages of the stress are trivial. 
On the other hand, the estimates for the standard deviations of the stress components, which follow from 
(39), are rather complicated and are omitted for conciseness. 

In the ideally plastic limit (m ^ 0), Eq. (45) may admit more than one solution for k varying with /. 
The solution of interest is the positive one (for / > 0) continuously connected to A: = 1 for m = 1, which 
solves 



{l-k) (efi'Vee) 
21 



k 



(49) 



The corresponding effective flow stress follows from (46) with m = 0. 

Dilute limit. We now evaluate the above expressions as / ^ 0. We begin by noting that from numerical 
solutions of (45) it can be inferred that k ^ m for both PS and SS loadings. 

When m > 0, the anisotropy ratio k remains strictly positive. In this case, Willot et al. (2008a) have 
shown that the effective modulus A appearing in expressions (42) has an expansion of the form (see expression 
(25) in that Ref.) 

A/A = l + Ci(A:)/ + C2(A;)/2 + ... (50) 

where Ci{k) = — (1 + V^), and C2{k) is assumed to have bounded first derivative and presumably depends 
on the specific loading. To first order in /, the strain variables (43) are then 

4'^ 1 . 1-k fi^^^V 1 s^^^ 1 

Se V2A:3/4^ ^2A:3/2^' J 2y^V2^' Se k^^^' ^ ^ 

independent of C2{k). Making use of these expressions, it is easy to verify that the zeroth-order term of the 
left-hand side of equation (45) is (1 — A:), which confirms that k = m -\- 0{f) as / 0. Finally, it is easy to 
verify that the expansion of expression (46) for ao, to first order in /, is given by (13) with (14). 

When m = 0, the anisotropy ratio /c ^ as / ^ 0. Numerical solutions of equation (49) show that 
k ^ hp^^. In this case, Willot et al. (2008a) have shown that the expansion (50) for the effective modulus 
A remains valid. To order /^/^, the strain variables (43) are then 

£e 1 + V263/4 (1 + V263/4)2 ' \^ ^ (1 + V2&3/4-)2 •' ' £^ 1 + ^263/4 

(52) 

The coefficient h is obtained by expanding equation (49) to order /^/^; the solution is 6 = 1/4, independent 
of C2(0) and in agreement with numerical results. Having determined the strain variables (52), it is easy to 
show that the expansion of expression (46) for a^^ to order /^/^, is given by (17) with ago = 1/2. Note that 
the result is independent of the second-order coefficient C2{k). 

A . 2. Stress formulation 

A dual version of the 'second-order' method follows from exactly analogous arguments. When specialized 
to porous materials with an isotropic matrix phase, the 'second-order' variational estimate for the effective 
stress potential is given by (Ponte Castaheda, 2002a; Idiart et al., 2006a) 

u{W) =stat |Sl (o^; ^^^] M^i) ) - (1 - f)v^^^ [a^^] M^^) ) | . (53) 

In this expression, ul is the effective potential of a porous LCC with a matrix characterized by a second- 
order, Taylor-type expansion of the nonlinear potential (2)2 about a reference stress tensor a^^^ given 
by 

4^\a;<TW,MW)=MW(<TW) + ^(aW):(a-aW) + i(o--a«):M«:(^ (54) 

OCT 2 

where M^-"^^ is a symmetric, fourth-order tensor (of compliances) of the form 

M(i) = ^E|| + i-E^, k = X/ii. (55) 
Here, A and /x are two shear moduli, and 

E|| = 2Tfo®TfLT' Ex=K-E||, (56) 
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are two orthogonal, fourth-order, projection tensors with principal axes 'aligned' with a^^^ 
In turn, the 'error' function v^'^^ in (53) is defined as 



^;(i)(^(i)M(i)) = stat{4^\a(i);,T(i),M(i))-«(i)(,T(i))} , (57) 

and is a measure of the degree of nonlinearity of u^^\ 

Given that the matrix potential (54) corresponds to a 'thermoelastic' material with a 'thermal strain' 
and a 'specific heat' given by 

^(1) = ^(^(i))_M«:a«, (58) 

ftd) = w(i)(d-W)-r,«:a«-i^W:M«:a« (59) 

the effective potential ul can be written as (Willis, 1981) 

ULiW) = ^a:M:W-\-rj:a-\-h, (60) 

where rj = r]^^\ h = {1 — f) h^^\ and M denotes the effective compliance tensor of a porous material with a 
purely elastic matrix with compliance tensor (55). It is useful to recall at this point that the first and second 
moments of the stress field cr 2, in the matrix of the porous LCC can be obtained from relations completely 
analogous to (31) (Idiart and Ponte Castaheda, 2007a). 

In addition, the reference stress a^^^ must also be specified in order to characterize the potential u^l^ . 
As already mentioned, an 'optimality' condition for this variable, if it exists, is not yet available, and an ad 
hoc prescription is therefore required. Following Idiart et al. (2006a) we identify &^^^ with the macroscopic 
deviatoric stress, i.e., 

^^^^=-^d. (61) 

Then, the stationarity conditions in (53) and (57) yield a system of nonlinear algebraic equations for the 
variables a^^\ A and /i. With (2)2, (55) and (61), the conditions resulting from the stationary operation in 
(57) read: 

(aJ^)-ae)/(3A)=^'(aW)(ai^)M'^)-^'(^e), (1/3^.) = (62) 



where cf^^\ = {?t/2)d'^^^ -^11,-1 •^^"'"^ represent the components of a^^^ that are 'parallel' and 'perpendicular' 

to the macroscopic stress o". The equivalent part of a^^^ is then (cf^^Y — + turn, the 

conditions resulting from the stationarity operation in (53) are 



.(1) _ , / 6 duL .(1) / 6 duL 

where a choice of roots has been made. Idiart et al. (2006a) Finally, making use of the relations (62)- (63), 
the estimate (53) can be written in the simplified form 

u{a) = (1 - /) [piaW) - V'(ae) (aj^^ - a^^)] , (64) 

where is the equivalent part of the average stress in the matrix of the LCC. 

Effective behavior. Differentiating expression (53) with respect to a yields (Idiart and Ponte Castaheda, 
2007a) 

e=g(a)=6^ + (l-/)7^i), (65) 
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where Sl = M<t + ?7 is the macroscopic strain in the LCC, and the (incompressible) tensor 7^-'^^ is given by 



7 



(1) - 



(i)_Mf^ 



L ) + 4 



3 A-i -/i-i 



(66) 



d'^u^^\cr) / dcrdcr is the tangent comphance tensor of the matrix phase evaluated at cr. 



Here, .v.^ 

Field statistics. Corresponding estimates for the first and second moments of the fields in each phase 
corresponding to (64) are given by (Idiart and Ponte Castaheda, 2007a) 



(s (g) e) 



(r) 



(r) 



f(0 



(<T a) 



(r) 



(67) 
(68) 



where s^"*^^ = du^^\a)/dcr^ in the porous phase s^^^ = 7^^^ = 0, and the subscript 5 in the tensor product 
denotes symmetrization. Thus, while the nonlinear estimates for the stress statistics in the matrix coincide 
with those quantities in the LCC, the estimates for the strain statistics exhibit certain correction terms 
which depend on the matrix nonlinearity through the tensor (67). 

Expressions for power-law materials under 'aligned^ loadings. When the macroscopic loading is 'aligned' the 
tensor M is given in terms of A and fi by the inverse of (40). Then, evaluating the derivatives of 2^ in (63), 
we arrive at the following relations: 



a5^V^e = l + yr|| + (l-aW/ae) , a^^^ /a, = f'f, 4^)^ = 1/(1-/), 
where the functions fii and f± are given by 



(69) 



(1-/) 



A-i 



1-/ 



(70) 



Again, the prime in these expressions denotes differentiation with respect to /c, and the ratio A/A is a 
function of k and /, and depends on the specific loading. Explicit expressions have been given in Willot et 
al. (2008a). Then, by taking the ratio between the two conditions (62) with ip given by (5)2, we obtain the 
relation 



In 



(1-lA) (a|fV^e)+(lA) 



In 



= 1 — n. 



(71) 



where ai^^ = ^ (<^|j^^)^ + (aj_^^)^. This relation, together with (69), constitutes a single nonlinear, algebraic 
equation for the anisotropy ratio k that must be solved numerically. 

Finally, the 'second-order' estimate (64) can be written as (6)2 with the effective fiow stress given by 



^ = (1-/)- 



1+n 



^ I -(l + n)| d^-^ 



(72) 



where and cr^^^ are given in terms of k by relations (69). The corresponding estimates for the field 

statistics follow from relations (67)- (68). The estimates for the standard deviations of the stress components 
in the matrix phase, as defined by (12), are given by 



5i?(l)(<T||,^) 

24 



(73) 



and the tensor 7*^^^ required to compute the strain statistics simphfies to 

7(^) = ^(A-^-Ar^)(aW-aW)^./ae, (74) 

where (2At)~^ = Ey -M^^^ = (3/2)?/;' (ag). However, the estimates for the strain statistics, which follow from 
(67) and (68), are rather complicated and are omitted for conciseness. 
In the ideally plastic limit (n oo), equation (71) simplifies further to 

^J'V^e = 1/(1 - k), (75) 
and the 'second-order' estimate (72) reduces to 

^=(a«/a.)"\ (76) 

(7f) \ / 



Dilute limit. We now evaluate the above expressions as / ^ 0. As in the strain version, numerical solutions 
of (71) indicate that A: ^ m for both PS and SS loadings. 

When m > 0, the anisotropy ratio k remains strictly positive and the eff"ective modulus A is given by 
(50) to second order in / (Willot et al., 2008a). To first order in /, the stress variables (69) are then 

J_ = ,,L^,./.^0(/n =5jL_/. (77) 

independent of C2{k). Making use of these expressions, it is easy to verify that the zeroth-order term of the 
left-hand side of equation (71) is (1 — l//c), which confirms that k = m-\- 0{f) as / ^ 0. Finally, it is easy 
to verify that the expansion of expression (72) for ao, to first order in /, is given by (13) with (14). 

When m = 0, the anisotropy ratio ^ as / ^ 0. As in the strain version, numerical solutions of 
equation (75) show that k ^ hp^^^ and the expansion (50) for the eff'ective modulus A remains valid. To 
order /^/^, the stress variables (69) are then 

The coefficient h is obtained by expanding equation (75) to order /^/^; the solution is 6 = 1/2^/^, in 
agreement with numerical results. Having determined the stress variables (78), it is easy to show that the 
expansion of expression (72) for a^^ to order /^/^, is given by (17) with aso = 3/2^/^. 

A. 3. Estimates for rigidly-reinforced composites 

Making use of the duality relations of Section 2.3, second-order estimates for rigidly-reinforced composites 
can be generated directly from the above estimates for porous materials. According to those relations, see 
Table 1, we have the following correspondences between the two sets of estimates: 1/m ^ m, {^o/ao)~^^'^ 
(ao/cro)^ SS (PS) loading PS (SS) loading, 1/k ^ k, strain (stress) version stress (strain) version, strain 
(stress) statistics stress (strain) statistics. 



B. Limit analysis bounds for periodic materials 

In the limiting case of a rigid-ideally plastic matrix, use can be made of the kinematic and static theorems 
of classical limit analysis to compute upper and lower bounds, respectively, for the effective ffow stress of 
the composite (Bouchitte and Suquet, 1991). The computation of such bounds amounts to evaluating 
the effective potentials (3) at (simple) trial fields 6(x) (or displacement u(x)) and cr(x), which in perfect 
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Figure 7: Trial fields used to obtain limit analysis bounds for porous materials, (a)-(c) Arrows denote displacement vectors, 
all of the same magnitude; (d)-(e) Arrows denote traction vectors associated with the stress field, all of the same magnitude. 
Dotted lines indicate surfaces of discontinuity. 
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(a) SS: < / < f (b) PS: < / < f 



Figure 8: Trial displacement fields used to obtain limit analysis bounds for rigidly-reinforced materials. Dotted lines indicate 
surfaces of discontinuity. 



plasticity include piecewise continuous fields with surfaces of discontinuity not only at material interfaces 
but also within the individual phases. Thus, this allows the use of piecewise uniform trial fields, which 
can lead to simple analytical bounds, as shown below. For consiceness, the procedure for computing these 
bounds is only briefly addressed here, and the interested reader is referred to Salengon (1983); Taliercio 
(1992); Maghous (1991). 

B.l. Porous materials 

Simple shear. Making use of the trial displacement and stress flelds shown in parts (a) and (d) of Fig. 7, 
the following bounds for 5o ^ire obtained: 

1 - 2727^/^2 < < 1 _ (2/V^)/V2. (79) 

Pure shear. Making use of the trial displacement and stress flelds shown in parts (b), (c) and (e) of Fig. 7, 
the following bounds for ao are obtained: 



VTT CTo 



i-J-f/\ 0</<f 

/f^ (80) 



Unlike in the case of simple shear loading, two difl'erent displacement trial fields have been used in this case, 
motivated by the FFT results shown in Table 2. The trial field shown in Fig. 7b yields the smaller bound 
for low to moderate porosities, while that of Fig. 7c yields the smaller bound for porosities closer to the 
close-packing limit. Both bounds coincide for / = 7r/6 ?^ 0.524. 

B.2. Rigidly-reinforced materials 

Figure 8 shows two displacement trial fields for SS and PS loadings. Both cases involve straight shear 
bands passing entirely through the matrix, from which it follows immediately that the resulting upper bound 
for 5q under both loading conditions is 

5o/^o < 1. (81) 
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This upper bound agrees exactly with the corresponding Reuss lower bound, and consequently it is the 
exact result. It should be emphasized, however, that in the PS case, the simple trial field of Fig. 8b is valid 
for / < 7r/8. For larger values of /, trial fields involving straight shear bands are no longer possible, and 
therefore the resulting upper bound on 5o will exhibit a finite reinforcement eff'ect. 
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